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Many phenomena of strongly correlated materials are encapsulated in the Fermi-Hubbard model 
whose thermodynamic properties can be computed from its grand canonical potential. In general, 
there is no closed form expression of the grand canonical potential for lattices of more than one 
spatial dimension, but solutions can be numerically approximated using cluster methods. To model 
long-range effects such as order parameters, a powerful method to compute the cluster’s Green’s 
function consists in finding its self-energy through a variational principle. This allows the possibility 
of studying various phase transitions at finite temperature in the Fermi-Hubbard model. However, 
a classical cluster solver quickly hits an exponential wall in the memory (or computation time) 
required to store the computation variables. Here it is shown theoretically that the cluster solver 
can be mapped to a subroutine on a quantum computer whose quantum memory usage scales 
linearly with the number of orbitals in the simulated cluster and the number of measurements 
scales quadratically. A quantum computer with a few tens of qubits could therefore simulate the 
thermodynamic properties of complex fermionic lattices inaccessible to classical supercomputers. 


I. INTRODUCTION 

The Fermi-Hubbard model (FHM) [1] is a central tool 
in the study of strongly correlated electrons in condensed 
matter physics [2]. It captures the simplest essence of 
the atomic structure of materials and the second quan¬ 
tization of the many-body interacting wavefunction and 
can be used to model phase transitions in Mott insula¬ 
tors, high-Tc superconductors EHH, heavy-fermion com¬ 
pounds [5], atoms in optical lattices, organic materi¬ 
als and many others. The exact solutions to the one¬ 
dimensional Hubbard model are known and well under¬ 
stood Hi] but the two- and three-dimensional models 
are known not to have general closed form solutions and 
are subject to important theoretical studies IMS]. An 
elegant approximation method valid for short-range in¬ 
teractions is cluster perturbation theory (CPT), where a 
lattice is divided into manageable identical clusters which 
are solved and then recomposed into a lattice through 
with perturbation theory [nun]. However, the method 
is not sufficient to systematically account for broken sym¬ 
metries in the FHM and has to be extended. In super¬ 
conductors and antiferromagnets, local interactions can 
have long-range effects and order parameters can appear 
in different regions of phase space. These effects can 
be taken into account in the Green’s function of a clus¬ 
ter by finding the stationary point of the lattice’s grand 
canonical potential when the self-energy of a cluster is 
taken as the variational parameter [15]. This self-energy 
functional theory (SFT) is a great computational tool to 
study the important macroscopic thermodynamic phases 
of the Hubbard model starting from its microscopic de¬ 
scription. In the context of the SFT, the CPT approxi¬ 
mation is generalized to what is known as the variational 
cluster approximation (VGA). 

However even simulating a small cluster with a hand¬ 
ful of electrons (or orbitals) is a difficult task for classical 
computers since the matrices involved in the computation 


scale exponentially in size with respect to the number of 
electronic orbitals. The quantity of information involved 
in the precise numerical treatment of large strongly cor¬ 
related electronic systems quickly reaches magnitudes 
where no reasonnable classical memory technology is suf¬ 
ficient to store it all. Therefore, being given access to a 
large controllable Hilbert space in a quantum computer 
offers the possibility of simulating electronic systems at 
the microscopic level with a greater complexity and ac¬ 
curacy than the ones accessible to classical computers 

US). 

This work is inspired from recently developed ap¬ 
proaches in quantum simulations such as the simulation 
of spin systems [nniHi, fermionic systems and quantum 
chemistry dSHm and boson sampling to extract vibronic 
spectra m- In general, it happens that the occupa¬ 
tion state of an electronic orbital can be efficiently repre¬ 
sented by one qubit on a quantum computer through the 
Jordan-Wigner transformation. The memory bottleneck 
in numerically representing the many-body wavefunction 
is overcome by making sure that it is never measured 
and stored on a classical memory at any point during 
the simulation. In the VGA, the quantities that need 
to be extracted from the wavefunction are the intra¬ 
cluster single-particle correlation functions whose num¬ 
ber scales quadratically with the number of orbitals in a 
given cluster. On the practical side, it is not yet known 
how the computing power of quantum processing devices 
will scale in the future, but machines with a fews tens 
or hundreds of qubits could already be very useful to run 
quantum subroutines as part of larger classical simulation 
algorithms. This proposed method could open a practi¬ 
cal way to model and engineer the electronic behavior 
of strongly-correlated materials with intricate crystalline 
structures in a unified and consistent manner. Further¬ 
more the underlying SFT is very general [23l [24] and 
not restricted to the class of FHMs. Similar schemes to 
simulate spin systems, the Bose-Hubbard model or more 
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exotic fields in lattice gauge theories Egni can likely 
be constructed in a similar fashion. 


This paper aims at at being self-containend by provid¬ 
ing all the concepts required to implement the solver on 
a general purpose quantum computer m It is struc¬ 
tured in the following manner. Section summarizes 
the variational cluster method used to compute proper¬ 
ties of the FHM. In subsection |II A[ a variational prin¬ 
ciple of the self-energy for the grand canonical poten¬ 
tial of the model is outlined such that it can account 


for possible long-range ordering effects. Subsection IIB 


formalizes the approximation where the Fermi-Hubbard 
lattice is divided in independent clusters linked with hop¬ 
ping terms. Section m introduces the detailed formal 
description of a cluster using the example of a 2D lat¬ 
tice with superconductivity starting in subsection [III A 


Subsection |IIIB| proceeds with reviewing the formalism 
to compute the Green’s function of the lattice from the 
independent clusters and subsection [III C| lists methods 
to compute observables of interest once the variational 
problem is solved. Section |IV| covers the computer in¬ 
tensive step where the eigenvalue problem of the clus¬ 
ter Hamiltonian must be solved at each iteration of the 
variational solver. Subsection IIV Al summarizes the so¬ 
lution method on a classical computer and a memory 
efficient quantum subroutine to introduced in subsection 


IV B The procedure to measure the Green’s function of 


the cluster is described in subsection IV C[ Appendix 
presents numerical results where the quantum procedure 
to compute a cluster’s Green’s function is shown to be 
equivalent to traditional solution methods. In appendix 
[B} details of the initial Gibbs state preparation are given 
for a specific algorithm. 


II. SOLVING THE FERMI-HUBBARD MODEL 
WITH THE VARIATIONAL CLUSTER 
APPROXIMATION 


The goal of this section is to introduce the impor¬ 
tant physical quantities of the main loop of the numer¬ 
ical variational solver used to extract properties of the 
FHM. Since the interesting observables typically corre¬ 
spond to the response of the system to external perturba¬ 
tions, the central object of study is the Green’s function 
which contains both the thermal and the dynamical prop¬ 
erties of the system. To compute the Green’s function, 
a variational principle on the grand canonical potential 
is derived from functional arguments. The Green’s func¬ 
tional variational problem is then mapped to a self-energy 
variational problem to account for possible spontaneous 
symmetry breaking from long-range ordering in a self- 
consistent manner. At last the lattice approximation is 
introduced to complete the description of the lattice vari¬ 
ational solver. 



Figure 1. Diagrammatic representation of the relation be¬ 
tween the single-particle Green’s function G, the bare Green’s 
function of the non-interacting lattice Got and the self-energy 

E. 


A. The grand canonical potential as a functional of 
the self-energy 

Variational solvers m are powerful tools to solve 
many-body problems in quantum mechanics. The FHM 
is an effective description of the microscopic physics of 
the electrons in a solid useful in calculating the proper¬ 
ties of Fermi liquids, Mott insulators, anti-ferromagnets 
[28], superconductors [29| and other metallic phases. The 
model describes a simple electronic band in a periodic 
lattice F where electrons are free to hop between orbitals 
(or sites) with kinetic energy t and interact via a simple 
two-body Goulomb term U. The standard form of the 
Fermi-Hubbard Hamiltonian is given by 

1-L t ^ ^ U ^ ^ ^ ^ (1) 

{i,j),a i i,o- 


where /i is the chemical potential that determines the oc¬ 
cupation of the band. The Cicr(ct^) are the fermionic an¬ 
nihilation (creation) operators and the number operators 
are riia- = Note that in the rest of this document, 

units are used such that t ^ 1 is assumed to be the ref¬ 
erence energy and inverse time. It is also assumed that 
h ^ 1 and ks 1. 


1. The Luttinger-Ward formalism 

The Green’s function G of the full system described by 
TL can be obtained exactly from the bare single-particle 
Green’s function of the non-interacting lattice (tight- 
binding) Got and the self-energy 5] = by 

solving the Dyson equation represented in figure 


G = Got + Got5:G. (2) 

When there is no interaction, the self-energy is zero and 
the tight-binding Green’s function for a given one-body 
hopping matrix t is 


— u; - t. (3) 

The model can be considered “solved” once the single¬ 
particle Green’s function G can be computed accurately 
for any interesting input coordinates (such as position 





















3 


G 



Figure 2. The Luttinger-Ward functional $ is the sum of all 
the two-body skeleton diagrams. The functional derivative 
with respect to G gives all the diagrams for the computation 
of Xl. In the case where U = 0 then $ [G] = 0. 


/ momentum, time / energy). A method to obtain the 
Green’s function consists in rewriting the Dyson equation 
as a variational principle on the grand canonical poten¬ 
tial of the system. To accomplish this task, it is useful 
to introduce the Luttinger-Ward functional HSl |30] of 
the Green’s function ^ [G] which generates all two-body 
skeleton diagrams (see figure and has the interesting 
property that its functional derivative with respect to G 
is simply 


6^ [G] 


= E. 


( 4 ) 


Furthermore, the functional form of ^ [G] depends only 
on the form of the interaction U and is independent of the 
one-body terms in 'H. In statistical mechanics, observ¬ 
ables are derived from a thermodynamic potential. For 
many-body systems, it is typically easier to let the total 
number of particles fluctuate and work with the grand 
canonical ensemble. The grand canonical potential of 
the full lattice can be defined from the Luttinger-Ward 
functional as a functional of G 


fit [G] = $ [G] -IV [(Go/ - G-i) G] +Tr In [-G], (5) 

such that the Dyson equation © can be recovered as the 
stationary point with respect to the variation of G: 


(5fit [G] 
(5G 


= S-Go/ + G-i = 0. 


(6) 


In Ref. [3T], Potthoff describes three types of approxi¬ 
mation stategy to solve this variational problem. A type 
I approximation would try to simplify the Euler equation 
from a heuristic argument but could suffer from thermo¬ 
dynamic inconsistencies. A type II approximation would 
correspond to computing the $ [G] functional only for 
a finite set of diagrams, but justifying the use of a par¬ 
ticular functional form over other possibilities is in itself 
not trivial. Finally, in a type III approximation, ther¬ 
modynamical consistency is preserved as well as the ex¬ 
act form of the Luttinger-Ward functional but the trial 
Green’s functions are chosen from a restricted domain 
where the self-energy is constrained. The VGA is a type 


III approximation. The main advantage of this type of 
scheme is that it allows for a systematic construction of 
increasingly accurate solutions to many-body problems 
with local interactions. In the case of the FHM, a good 
scheme to systematically approximate the self-energy is 
to consider a reference lattice of isolated clusters F' with 
the same local interaction term U as the lattice F and 
pick 5] from the exact solution of the reference lattice. 
This method allows for the construction of solutions to 
the FHM that are very accurate except for long range 
correlations that exceed the dimensions of the clusters. 
The main advantage of this scheme is that the solutions 
are guaranteed to become asymptotically exact as the 
size of the cluster reaches the size of the original lattice. 
The next step consists in rewritting the grand canonical 
potential Qt as a functional of the lattice self-energy E 
instead of the Green’s function G. 


2. Self-energy funetional theory 

The variational principle of the self-energy of a clus¬ 
ter [32] intends to account for solutions of the Hub¬ 
bard model with spontaneous symmetry breaking caused 
by long-range interactions. The grand canonical poten¬ 
tial Vtt [G] can be rewritten as a functional of the self¬ 
energy rtt [^] by applying the Legendre transformation 

G [S] = (Go/ - S)“^ such that 

fit [G] = $ [G] - IV [(Go/ - G-i) G] + IVln [-G] 

= $ [G] - Tr [SG] + Tr In [-G] 

'-V-" 

A[S] 

= A[S]-IVln[-Go/ + S] 

= fit[S]. 

( 7 ) 

Let’s then notice that [^] is still exact and now only 
depends on the self-energy E and the non-interacting 
Green’s function Got- The Legendre transformed 
Luttinger-Ward functional A [E] has the nice property 


M [E] 


( 8 ) 


which is used to recover the Dyson equation of the system 
and the variational principle depending on the self-energy 


^^^ = (Go-/-S) '-G = 0. (9) 

Solutions to the FHM can be found by varying the self¬ 
energy until a physical value of the Green’s function is 
found and the Dyson equation is satisfied. However, since 
this is in general a saddle-point problem, the optimal 
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Figure 3. The essence of the VGA method is to remove the 
one-body links (denoted t) between small clusters (contained 
in V) from the lattice F and consider only the reference lat¬ 
tice F' whose Hamiltonian Ti' is block diagonal in the Wannier 
basis and easier to solve than the complete problem 1-L. The 
reference system generates a manifold of trial self-energies E' 
parametrized by single-particle parameters t^ The self-energy 
functional can be evaluated exactly on this manifold as the in¬ 
teraction part of the Hamiltonian (the Us) is left unchanged.. 
The solution become asymptotically exact as the clusters are 
made to include more sites. 


point cannot be interpreted as an upper bound to the 
exact energy (as in the Ritz variational method) but as 
the most “physical” approximation of the grand canonical 
potential allowed by a given parametrization of the self¬ 
energy. Computing the exact single-particle self-energy 
for a large lattice and storing the result are tasks be¬ 
yond the capabilities of classical computers. The idea of 
cluster methods used to approximate the solution of the 
full lattice F is to divide it into a reference lattice F' of 
clusters of a small number (i.e. computer tractable) of 
sites, solve a cluster exactly and use perturbation theory 
to approximate the properties of the full lattice. 


B. Variational cluster approximation 

Large lattices with millions of orbitals are impossible to 
simulate exactly on classical computers since the memory 
required to store for the associated state vectors scales 
exponentially in cluster size. A method to mitigate this 
problem makes use of the translation invariance of the 
lattice. It consists in breaking down the lattice in several 
independent clusters and making use of the universality 
of the Luttinger-Ward functional to recast the variational 
equation on a cluster-restricted domain of the self¬ 
energy. The exact solutions are recovered when the size 
of the cluster is equal to the size of the original lattice 

m- 

Good and thorough introductions to the VGA method 
can be found in [HISll. In the restricted Hilbert space 
of a cluster, the goal is to variationally find a self-energy 
E' such that it is most physical (by satisfying the VGA 
version of the Dyson equation) and minimizes the free en¬ 
ergy. As hinted at the end of subsection |II A| and shown 
in figure the VGA approximation consists in subdi¬ 


viding a full lattice F into a reference lattice of identical 
clusters F' and solving the reference model exactly in 
order to obtain its self-energy E'. In this context, the 
Green’s function of a cluster is a frequency dependent 
matrix given by 

G'-^(a;) =a;-t'-E'(cj) (10) 

The Legendre transformed Luttinger-Ward functional A 
only depends on the interaction part of the Hamiltonian. 
Since by definition the interaction part of the Hamilto¬ 
nian is the same for the full system and the reference 
system, the identity A [E)'] = A [E)] must hold. Let’s 
note that this scheme would not work directly in the 
case of the extended FHM (where there is intersite in¬ 
teraction), since a reference system of independent clus¬ 
ters cannot be found by simply removing one-body links 
of the Hamiltonian [23]. As in equation 0, the grand 
canonical potential of the reference system is given by 


n' = Qt' p'] = A p] - Tr In [-G'] , (11) 

where G' is the Green’s function of the reference sys¬ 
tem. When they are both evaluated at the self-energy 
of the reference system, the difference between the grand 
canonical potential of the full lattice and the reference 
system is 


Qt P'] =n' + Tt In [-G'] - Tr In [-G] . (12) 

This relation is exact, the only approximation of the VGA 
is in the restriction of the domain of the self-energy. It 
can be further simplified as the VGA is built within SFT 
as a well-defined variational extension to the CPT. The 
full lattice Green’s function G [E)] is equal to the CPT 
Green’s function if its self-energy is restricted to the do¬ 
main of the reference system. As in figure it is useful 
to define V = t—t' as a perturbation, where t contains 
all the one-body terms of the full lattice F and t' repre¬ 
sents all the one-body terms of the lattice of clusters F'. 
As a result of strong-coupling perturbation theory, the 
CPT Green’s function is given by 

G[S'] = Gept = (G'-i-V)-\ (13) 

With some algebra, equation (|12|) can be written as 


rit [S'] = fl' - Tr In [1 - VG'j. (14) 

The functional is exact as no classes of diagrams have 
been explicitly excluded. At the saddle-point, it repre¬ 
sents the quantity which is physically the closest to the 
physical grand canonical potential of the full lattice when 
the self-energy is computed on the reference lattice. The 
effect of single-particle correlations and intra-cluster two- 
particle correlations is treated non-perturbatively but the 
inter-cluster two-particle effects are neglected in the one- 
particle spectrum. Even if only a small cluster is exactly 
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solved, the self-energy variational principle can be 
used to study the properties of the infinite system like the 
various order parameters in a thermodynamically consis¬ 
tent framework. Since the VGA is a well defined general¬ 
ization of the CPT, it also shares similar characteristics. 
It is exact in the limit y ^ 0 where the self-energy dis¬ 
appears to yield the tight binding model. It is also exact 
in the strong-coupling limit ^ ^ 0, where all sites are 
effectively decoupled. The method is easy to generalize 
to non-homogenous lattices. The next section introduces 
the details of the objects required to compute (14) and 
find its stationary point as well as some observable that 
can then be calculated. 


III. EXAMPLE ON A SQUARE LATTICE WITH 
SUPERCONDUCTIVITY 

In this section the self-energy variational approach is 
used to model superconductivity in a Fermi-Hubbard lat¬ 
tice. A more general formulation of possible orders could 
be made (for arbitrary ordering potentials and cluster 
graph), but the goal of this section is only to introduce 
the types of formal elements required to describe a clus¬ 
ter. Other types of order parameters can be found in 
the literature [5]. First the different terms in the Hamil¬ 
tonian of the cluster are explained. Then the detailed 
formalism of the VGA is given through the example of 
a square lattice with superconductivity. Finally, various 
quantities involved in the computation of useful observ¬ 
ables are listed. 


A. Hamiltonian of a cluster 

Each cluster includes only a small portion of the terms 
of the original lattice and variational terms must also be 
included to account for possible long-range order. For 
convenience, let’s assume that F is a square lattice with 
constant spacing a. It is broken down into Nc clusters 
each with Lc orbitals (“sites”) with two electrons each 
(spin up t and spin down ^). The Hamiltonian of each 
cluster is given by 


H' = HfH + ^local + ^s-pair + ^d^ 2 _y 2 + ^AF, (15) 

where the Fermi-Hubbard terms remaining in F' are given 
by 


^local — M ^ ^ ^ 


(17) 


It can be seen that at the stationary point 
electronic occupation expectation value is 


/V ^ ^ 

(n) = Tr G =-;— = 

a/i 


dflt dflt dt' 
dfi ^ dt' dji 


0, the 

(18) 


where the two methods converge to the same average 
occupation at the stationary point. Keeping the chemical 
potential fixed in the cluster Hamiltonian would break 
this condition. 

The spontaneous transitions of the FHM can be stud¬ 
ied by introducing artificial symmetry breaking terms to 
the cluster Hamiltonian and treating them as variational 
variable. The choice of these terms is somewhat arbi¬ 
trary and is usually justified by the physics of the system 
studied. For example in the FHM, it is often interesting 
to study the competition between superconducting order 
parameters with different symmetries and the antiferro¬ 
magnetic ordering. A variational singlet pairing term is 
introduced as 


^s-pair — A 5 (19) 

i 

while a dj. 2 _y 2 singlet pairing takes the form m 

^d^2_y2 = dij -h CjiCi^^ , (20) 


where R are the vector positions of the sites in the cluster 
and 


dij 



if Hi — Hj = ±aex 
if Hi — Hj = ±aey 
otherwise. 


( 21 ) 


The variational Neel antiferromagnetic Weiss field takes 
the form 


HAF = M'^e*Q-*^‘(nit-n4), (22) 

i 

where Q = (tt, tt) is the antiferromagnetic wavevector. 

The small parameter in the approximation is 
which means that increasing the size of the cluster also 
increases the accuracy of the simulation. 


^FH = E - U (16) 

{i,j),o- i 

which is the same as 0 without the chemical potential 
term. The chemical potential must be kept as a varia¬ 
tional term to enforce the thermodynamic consistency of 
the electronic occupation value 


B. The superlattice of clusters 

The relation between the original lattice and the lat¬ 
tice of cluster is given in more details along with useful 
notations. The main objects of interest for the quantum 
subroutine are introduced in this subsection. 
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(26) 


which is done efficiently on a classical computer once 
fit [t'] can be evaluated for a given set of parameters (for 
example, by a Newton-Raphson method). In the case of 
a lattice problem, the grand potential functional takes 
the following form 


Figure 4. Reduced Brillouin zone of the reciprocal lattice. 
The quasi-momentum vector k belong to the reciprocal lattice 
of r while the K component belongs to the reciprocal lattice 
of a single cluster. The k vector belongs to the reciprocal 
superlattice (hopping between clusters). 


1. The superlattice in reciprocal space 

To make the procedure clear and concrete, let’s work 
on the example of the superconducting order parameter 
on a 2D lattice. For a good explanation of quantum clus¬ 
ter theories and the details for computations on clusters 
of arbitrary size see OEll. A square lattice with 8 or¬ 
bitals per cluster is required to study s-wave and d-wave 
superconductivity in the FHM. Let’s take a lattice F with 
N sites and divide it in clusters of Lc = 2 x 2 = 4 sites, 
then the number of clusters is simply Nc = ^. Let’s 
label these 4 sites as 1, 2, 3 and 4. When the full 
lattice is Fourier transformed, the first Brillouin zone in 
quasi-momentum space is given by 




I-V 


(k) G' {z)_ , 


(27) 

where V ^k^ and G' (z) both depend on the chosen vari¬ 
ational parameters (the hat notation is explained below, 
it refers to the Nambu space). The contour integral dz 
can be done as a real line integral, as a Matsubara sum or 
as an efficient summation based on the continued fraction 
expansion of the Fermi function |35| . 


3. The eigenvalue problem 

In order to evaluate the energy-dependent Green’s 
function G' (z), the eigenvalue problem for one cluster 


n' 10,) = En |0n) (28) 


= m,,y = 0, (23) 

and the reciprocal superlattice is given by 


must be solved for different parameters until the station¬ 
ary point is reached. For 2Lc orbitals , the eigenvalue 
problem of the Hamiltonian can be solved in the occupa¬ 
tion eigenbasis defined by 


k 


xjy 


/ y 

Na ’ 


3.x/y 


AT, - 1 . (24) 


2. The saddle-point problem 


The observable properties of the Hamiltonian § can 
be computed from the CPT formula by finding vari¬ 
ational parameters (/i',A', M', etc.) that generate 

D for which the Dyson equation is stationary. In 
practice, this condition is reformulated explicitly over the 
variational parameters as 



(25) 


In the superconducting Fermi-Hubbard example, this 
would correspond to solving the saddle-point problem 


y-j / I \ Tlit y-y / ^ \ 

riL.i) = |Vac), 

i=l i=l 

(29) 

where |Vac) is the many-body vacuum. The dimension 
of this Hilbert space is 4^^ which means that storing the 
matrices of the calculation scales prohibitively with clus¬ 
ter size on a classical computer. Let’s note that spatial 
symmetries that commute with the cluster Hamiltonian 
TL' can be used to reduce the memory requirement of the 
computation HU- In all cases, it is useful to introduce 
the Nambu (singlet particle-hole) space notation. This 
notation is especially useful when considering quantum 
mechanical problems where an order parameter can ap¬ 
pear from broken gauge symmetries. In this space, field 

operators are replaced by a vector C 4 j such 
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that the energy-dependent Green’s function of a cluster 
can be represented in the form 


G'M = 


G' (w) F' (w) 

F't - G' ( - cj) 


(30) 


where the elements (w) = are the com- 

ponents of the single-particle Green’s function and 
Flj (uj) = {ci^Cji)^ are the components of the anoma¬ 
lous Green’s function. The (.. .)^notation corresponds 
to the frequency-dependent correlation function (i.e. the 
Fourier transformed two-point time correlation function). 
In the 4-site example, these matrices would have the form 


.t 


G' M = 


C2^c\^ 


C 2 tCit/ 


^2t4t 


t 


i^^C4tclt)^ ^C4t4t)^ (c4tT)^) 


and 


F'(w) = 


(ClfCu)^ (CltC24,)^ (CltC3t)„ (CltC44,)^ 
(C2tcu)^ (c2tC2t)^ (c2tC3i)^ {c2\Cii) ^ 
(C3tcu)„ (c3tC2t)t^ (c3tC3^)„ (c3tC4t)(, 
\ (C4tCu)^ (C4tC24,)^^ (C4tC3t)„ {ca^Ca^)^ 


(31) 


(32) 


Methods to evaluate G' (cj) on classical and quantum 
computer are given in section [IV| 

Let’s notice that in the Lc = 2 x 2 cluster, 32 different 
correlation functions have to be evaluated. In the general 
case, the number of correlation functions simply scales 
as 4-1/^, which is much smaller than the exponential 
scaling required for storing the full density matrix. See 
subsectior jlV A| and subsection |IV B| for the procedure to 
obtain these Green’s functions. 

At this point, the CPT potential in the reciprocal su¬ 
perlattice basis can also be defined as 


V (k) s t (k) - t', 


(33) 


where t contains all the one-body terms of the bare 

lattice r (i.e. no interaction terms) of the Hamiltonian 
0 . For the example of the square lattice, this gives 


where 


(34) 



_ 


(k). 



M e (fcj;) e (fcy) 0 ^ 


e K 


y 0 e* (kyj €* -n j 


(35) 


r 


and the dispersion relation for the square lattice is 


(fc) =-t(l + , 


(36) 


The t' term in equation (33) contains all one-body terms 


of a cluster (15), including the variational terms. In the 
example, 


t' = 


B C 
C D 


(37) 


where 


B = 


and 


D 


-/i' + M' 

—t 

—t 


0 

\ 

—t 

-li' - M' 

0 


—t 


—t 

0 

-li' - 

M' 

—t 


0 

—t 

—t 

-m' 


'/ 






(38) 

li' FM' 

t 

t 

0 

\ 


t 

t 

li' -M' 

0 li' 

0 

-M' 

t 

t 


(39) 

0 

t 

t 

li' FM' 

/ 



The pairing part is given by 


A' A' -A' 0 \ 


C = 


0 -A' 


a;; A' 

0 

0 -A' A' A' / 


(40) 




















The lattice-perturbed Greenes function 


Once the saddle point = 


m'* \ 

A' 

— ^ 


of equation (25) 


\ MU 

and V ^k, t* 


eval- 


is found, the function G' (cj,!* 

uated and the lattice-perturbed Green’s function can be 
computed. From here the dimensionality of the matri¬ 
ces involved in the calculations scales only as the square 
of the number of orbitals and can be performed easily 
on a classical computer. The lattice-perturbed Green’s 
function can be calculated to first order as 


g(k,Lo) = (G'-iH-v(k)) ' 


At this point the problem is solved and many observable 
quantities can be computed efficiently |28| . 


C. Calculation of observables 

Based on [36] , this subsection contains examples of ob¬ 
servables useful in explaining the result of experiments 
and landmark properties of the FHM. 

The average particle density is 




i=l 


and must agree with the value given by (18). The chem¬ 


ical potential /i can be scanned until a desired value of 
n is found. For superconducting problem in the FHM, 
it is useful to fix the chemical potential /r such that the 
lattice is maintained at quarter filling n = 0.25. The 
superconducting gap is given by 


A = UU = • (43) 

To recover the Green’s functions of the full lattice F, 
the “clustering” (which is a unitary transformation) is un¬ 
done and, taking into account the artificial translational 
symmetry breaking of the lattice, the single-particle and 
anomalous GPT Green’s functions are recovered in the 
lattice reciprocal space 


0cpt(k,u;) = 


o-ik-(ri-rj) 




( 44 ) 


From these quantities, the single-particle quasiparticle 
spectrum and the Bogoliubov quasiparticle spectrum can 
be evaluated as 


^(k,a;) = -piin^^o+Im^?cpt(k,w + ir/) 

F{k,Lo) = -piin^^o+Im7'cpt(k,w + i7?), 
from which the density of states is found to be 

AH = • 


(45) 


(46) 


The Fermion momentum distribution and the con¬ 
densation amplitude momentum distribution are respec¬ 
tively given by 


(47) 


For the case of a lattice with superconductivity, an inter¬ 
esting observable is the pair coherence length in real and 
reciprocal space given by 


( 41 ) 

A(k) 

~ 2772 ^cpt (k, z) 

X ic- 

F(k) 

~ (k, z) 


r = 


E,r^|A(r)|y _ Ei.|VkA(k)|^ 
Er|A(r)|^ Ek|A(k)|" 


(48) 


Depending on the problem, more observable can be 
computed with similar methods. Note also that the con¬ 
tour integrals map to the following form in the real time 
domain 


/ ^ 


(49) 


in the case where the retarded part of the Green’s func¬ 
tion is used to compute the integral. The Fermi function 
has the usual form f (uj) = -The self-energy 


1+e T 


variational approach has been outlined and the method 
which starts with a Hubbard-like description of the mi¬ 
croscopic details of a given solid and compute its ther¬ 
modynamic properties in a systematic way is complete. 
The next section reviews how the eigenvalue problem (28) 
is typically solved on classical computers and introduces 
the quantum subroutine. 


IV. 


SOLVING THE EIGENVALUE PROBLEM 
ON A QUANTUM COMPUTER 


Solving the eigenvalue problem (28) for a large num¬ 


ber of electrons is exponentially costly in memory as the 
number of orbitals increases. This section is divided the 
following way. First the classical eigenvalue solver for the 
Green’s function is described. Then the Jordan-Wigner 
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Number of orbitals 

Memory required 

3 

1 KB 

8 

1 MB 

13 

1 GB 

18 

1 TB 

23 

1 PB 


tude matrices then take the form 


(et) 


f Qlr 


Q = 


Qill 

Wlr 


(50) 


Table L Order of magnitude estimation of the classical mem¬ 
ory required to store the full finite temperature density matrix 
of a cluster with a given number of irreducible orbitals for a 
general cluster Hamiltonian. It is assumed that each matrix 
element is stored as a complex double-precision number (16 
bytes/element) and no optimization is used. 


transformation is used to map the cluster Hamiltonian to 
a quantum register. A method to generate initial Gibbs 
states in a quantum computer in reviewed and finally 
a procedure to extract the Green’s function out of the 
Gibbs state is explained. The full quantum procedure is 
shown to be efficient in quantum memory resources. 


A. The method on a classical computer 






/ 


and can be recast as Q' = QVft at non-zero tem¬ 
perature. It is also useful to define and compute 
g {uj) = It an be noted that Q is a 2Lc x 16^^ 

matrix which scales exponentially in memory with 
the size of the system being studied. 


5. Then compute (30) as G' (cj) = Q'g(co’)Q'^. This 
is the most time-consuming step on a classical com¬ 
puter, especially at non-zero temperature. 


6. The grand potential functional (27) and the lattice- 
perturbed Green’s function (41) can then be evalu¬ 
ated to respectively solve the saddle-point problem 
and compute observables. 


The resource intensive part of the numerical varia¬ 
tion solver is the computation of the energy-dependent 
Green’s function of the cluster G' (cj,t). On a classical 
computer, the memory used to store the description of 
the state of the system scales exponentially in system 
size. 


Typically, the Hamiltonian (15) is encoded in the oc¬ 
cupation basis (29) and the Schrodinger equation (28) 
is solved explicitly using an appropriate numerical diag- 
onalization method. As shown in table [Ij the memory 
usage scales exponentially with system size and diago- 
nalization typically scales as O (L^) in the number of 
arithmetic operation required. When successful, a set of 
eigenvalues {En} and associated eigenstates {|0n)} are 
obtained. If the cluster has Lc sites with 2 electrons each 
(spin up and down), then there are 4^^ eigenstates. The 
rest of the procedure is the following: 


1. Write ujmn — E^i E^yi. 

2. Write the occupation probabilities Pmn = 
e ^ ^ + e ^ ^ ^ Note that P = T~^ is the inverse 
temperature and Z = Tre“^^ is the partition 
function. 


3. Define the electron-like and hole-like amplitude 

Qimn ~ I Qt 1^^) Qimn ~ (0m I |0n) • 

4. Vectorize the m, n —r indices to obtain the ma¬ 
trices Ers = ^ fgUjf and Il'j^g — ^Ts Pr . The ampli- 


B. The method on a quantum computer 


Computing the Green’s function of the cluster G' (cj, t) 
on a quantum computer is possible in a hybrid analog- 
digital simulator. The first step generates a Gibbs state 
PGibbs (T) with some temperature T (or p = ^) measured 
on the digital register and the second step measures the 
correlation function of the cluster on an analog channel. 
The Jordan-Wigner transformation is used to map the 
Fermi-Hubbard Hamiltonian to a quantum register. The 
general procedure is the following: 


1. Map the cluster Hamiltonian (15) to a qubit system 
with the Jordan- Wigner transformation. 


2. Evaluate the two-point correlation functions 

for many different times for at least a full Hamil¬ 
tonian cycle (at zero temperature) or until corre¬ 
lations flatten out. Fourier transform to obtain 
the frequency-dependent correlation functions. The 
Hamiltonian is evolved in time using Trotter steps. 
Note that in the Jordan-Wigner basis, O (Lc) gates 
are needed at each time step. The full density ma¬ 
trix does not need to be measured, only O (L^) cor¬ 
relation functions need to be evaluated. 

3. Again, the grand potential functional (21^ and the 
lattice-perturbed Green’s function 0 can then be 
evaluated efficiently on a classical computer (simple 
linear algebra on small 2Lc x 2Lc matrices) to re¬ 
spectively solve the saddle-point problem and com¬ 
pute observable. 


























10 




Figure 5. Cir cuit to simulate the time-de pen dent correlation 
function (66) of the cluster Hamiltonian (15). The first part 
meant to generate a Gibbs state is taken from ED- Register R 
is used in the modified phase estimation scheme to prepare a 
rectangular state between the bath and the system contained 
in register Q. When the bath is traced out the system channel 
is left in a Gibbs state from which the different correlation 
functions can be read from the one-qubit register P. The 
size of register Q depends on the number of orbitals in the 
simulated cluster (typically n = 2Lc) and the bath size (which 
can be some constant factor larger than the system register). 
Register R is used as a digital component and q is therefore 
be the size required for the desired floating point accuracy on 
reading s^. Note that the numbers in the controlled gates of 
register R denote the index of the qubit which is acting as the 
control. 


mapping that requires no oracle black box for PL'. The 
Hamiltonian (15) is broken into M non-commuting parts 
such that 


M 

h' = J2k- 


(51) 


Each time-step Ar evolution of the cluster Hamiltonian 
m can be simulated with ut Trotter-Suzuki steps 


—il-L' At 


f M 

n< 

Vi=l 


H-L'-At 


E 

l<j 


[%[,%',] Ar 2 


2tit 




(52) 

The size of those time-steps set the upper bound in the 
simulated energy spectrum which scales as oomax oc 
while the lowest energy scales at the inverse of the total 
simulation time. 

The creation and annihilation operators of the Hamil¬ 
tonian can be mapped to the quantum computational ba¬ 
sis using a Jordan-Wigner transformation [42]. If there 
are 2Lc electrons, then the Jordan-Wigner m trans¬ 
formed creation operators are given by 


A full quantum circuit to measure G' (ci;,!) is shown 
in figure The specific algorithm [37] to create a Gibbs 
state was chosen mostly for aesthetic reasons. It appears 
to be the only Gibbs state generation method that pro¬ 
vides bounds on all parameters of the algorithm and that 
can be written in a circuit model. For completeness the 
main results of m are summarized and commented in 
appendix!^ There is no reason to believe that other sam¬ 
pling methods [38ti4Q] would not work also. A variational 
eigensolver m or an adiabatic quantum algorithm m 
could hypothetically be used to supply the initial ground 
state in the case of a simulation at zero temperature. 

Equation (28) does not need to be solved explicitly on 


a quantum computer, only a few correlation functions of 
interest need to be computed, this is explained in details 
in subsection |IV C[ The controlled evolution gates shown 
in figure assume that the Hamiltonian of the cluster 
can be mapped to a Hamiltonian in the quantum com¬ 
puter Hilbert space. Here is the procedure to make the 

I 


c\^ = I® 2 Le-i ^ 0 -+o erf*-! 

gt - 

In this notation, 


(53) 


I (7+ 0 






1 

(J 

cr G CT^ 


k-1 


A: = 0 
k = l, 
k>l 


(54) 


also <7+ = ^ and (Jz = 2(Jn — I, where 

(Jn = cr+cr_. The relations = —(Jz(J+ and 

(JzCF- = cr_ = —G-Gz can also be used. Note that 
the Jordan-Wigner transformation is independent of the 
Hamiltonian of the system and the dimensionality of the 
system. 


In the Pauli basis of the quantum computer, the terms of the cluster Hamiltonian (15) transform to 


^ + ^jcrGcr 




) 


U ni-fUii 

{4Ai + ciiCit) 


-t 0 Tl^ (i,j) + {i,j) 01®-'"=) 

-li' ( 1 ®-^= 0 Tl, (i) + Tl, (i) 0 1 ®-^=) 
U{Tl^ (i)0Ti, (i)) 

{i,i). 


( 55 ) 
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The strings of Pauli matrices are defined as 

{i,j) = (g) (8) (g) cr_ + cr_ (g) (g) cr+) (g) (56) 


where i > j between 1 and Lc and 

Tl, (i) = 1®-^=“* O cr„ ® 1®*“^ (57) 

Since and (i) conserve total spin in the Pauli basis, they are also number conserving in the occupation 

basis. For pairing terms it is also useful to define 


(i,j) = 1®-^=-^ o (cr+ O (jfLc-i+i-i (g) cr+ + cr_ (g) fjfLc-i+j-i ^ ^ 1®*“!. (58) 


In this case, i and j can be anything between 1 and Lc. The terms of D^c (bj) do not conserve total spin in the Pauli 
basis as they do not conserve the total number of particles in the occupation basis. 



Po Pt 


cuit is shown in figure It is a variation on DQCl 
(deterministic quantum computation with one quantum 
bit) 06] and phase estimation. 

A thermal density matrix of the simulated system must 
first be prepared in register S 


Figure 6. Circuit to measure the correlation function (66) 
from an input Gibbs state. Register S initially contains 
a given Gibbs state at inverse temperature /3 and register 
P is a single qubit initialized in the zero state. P is put 
in a state superposition by applying a Hadamard gate M 
and then used to apply the controlled evolution sequence 
(t) = Ul (r) a^Us (r) with Us (r) = ^ to the 

system channel. Finally the state superposition is reversed 
by a last Hadamard gate and the measurement in repeated to 
obtain the probability P(At), which returns information on 
the cluster Green’s function dSOl. 


In cases where the number of electrons in conserved 
in the cluster Hamiltonian (with superconductivity, the 
anomalous pairing terms break this symmetry), it is pos¬ 
sible use a Bravyi-Kitaev transformation m for an im¬ 
provement in the quantum memory usage of the algo¬ 
rithm (O (In Lc)). The mapping of PL' to the quantum 
computer is known and a method to generate Gibbs state 
has been chosen, the correlation functions can be mea¬ 
sured. 


Po = PGibbs W) G |0) (0|, (59) 


where 


PGibbs (/d) — ^ ^ ^ < 




|0m) (071 


(60) 


is a Gibbs state at some given temperature. It is to be 
expected that preparing a low temperature Gibbs state 
(large 0) is hard in general |47|, while high temperature 
Gibbs states 0^0 are simply fully mixed states which 
are easier to prepare. 

A sequence of controlled gates and controlled Hamil¬ 
tonian evolution follows the application of a Hadamard 
gate on register P. The unitary evolution generated by 
the cluster Hamiltonian (15) is defined as 


Us{t) = 


= {(t>m\ ■ 


(61) 


For convenience of notation (as seen in figure |^ , it is 
useful to introduce the set of gates 


C. Measuring the correlation function 

In this section an analog circuit is used to measure the 
correlation functions of a cluster Hamiltonian at some 
temperature T using a variation of the phase estimation 
algorithm is explained [44].. The Nambu single-particle 
Green’s function of the cluster G' (cj,t) can then be re¬ 
covered from the correlation function. The quantum cir- 


0^1/ (r) = ul (r) a^Us (r) (62) 

that define the application of a self-adjoint operator 
on the system (detailed below), followed by forward time 
evolution, then the application of another and finally 
a reverse time evolution. When applied to a Gibbs state 
in a phase-estimation circuit, the state of the computer 
at time r is described by 



























12 


Pt — 4 (pGibbs H“ PGibbs^]i,zy ^fiu ) PGibbs H“ ^/xz^ ) PGibbs^ji,zy )) ^ |0) (0| 

H“ 4 (pGibbs PGihhsOjxu ('^) ~l~ ^/xzx ('T) PGibbs ^fiu ij ') PGibbs^ji,zy ij ')) ^ 1^) I 

H“ 4 (pGibbs H“ PGibbs^]i,zy ('^) ^/xzx ('^) PGibbs ^/xzx ('^) PGibbs^ji,zy ('^)) ^ 1^) (^1 

H“ 4 (pGibbs PGihhsO^^y (t) (t ) />Gibbs H“ ^/xzx ) PGibbs^ji,zy ('^)) ^ I I • 


It can be seen that pr contains the information of the 
correlation function (cr^ (r) cr^y (0)), which can be mea¬ 
sured by evaluating the probability P^^y (Ai = 0(1) ,r) 
of measuring zero (one) in register P (and then Fourier 
transformed to obtain {a^a^)^). Formally, the interest¬ 
ing correlation functions that need to be extracted have 
the textbook form m 


Cf,^{T) = (cr^ (r) cr^ ( 0 )) 

= Tr [pGibbsOy (t) + (t) pGibbs] (64) 


SO low that the input Gibbs state is effectively is a pure 
(non-degenerate) ground state, it is plausible that adding 
qubits to register P would speed-up the measurement of 
the P^iy’s in the traditional sense of phase estimation [iH] . 
The retarded Green’s function can be computed numer¬ 
ically as 

(j) = -i9 (r) (t) (67) 

where 0 (r) is the Heaviside function. It can be Fourier 
transformed to get the Green’s function in the frequency 
domain 


_ p-iT{Em-En) Amn 

where A™” = e 

Note that these functions always outputs a real num¬ 
ber. If the controlled operation c — O^y (r) is applied 
for a time r > 0, the phase estimation algorithm yields 
the following probability for the two different outcomes 
M = 0 and = 1 


{M = 0,t) = i (1 + (t)) 


P^AM = 1,t) = i (1 - (r)). 


(65) 


Then from measuring the probability trajectory, the 
functions (64) can be recovered as 


/ oo 

dr (r). (68) 

-OO 

The spectral function can be obtained from the retarded 
Green’s function as 


Afiu (^) — ^ (cj) - (cj)) 

= -nm{G^^(w)}. 


(69) 


Since creation and annihilation operators are not Her- 
mitian, they cannot be used as and directly. A 
trick consists in using a linear combination of the oper¬ 
ators. For each electron orbital, the Hermitian Xi^- and 
Yia- operators are defined from (53) such that 


C^y (r) = 2 {P^y (M = 0, r) - P^y (M = I, r)). (66) 

As in DQGI [45], in general it is not useful to use mul¬ 
tiple ancillary qubits and an inverse Fourier transform 
to extract multiple bits of the probabilities P^y at each 
measurement shot since the input pGibbs is a state mix¬ 
ture. In the case where the simulated temperature is 


A^icr - Qcr T Grr 




(70) 


Note that T^cr 


transformation 




where Zjrr = cL 


\. The elements of (30) can be computed from the inverse 


^ Uia ('r)c]^, ( 0 )\ ^ 
(4a ('t) Cia' ( 0 )) 

(Cia ("t) Cja' ( 0 )) 

V(4a(^)ck ( 0 ))/ 




I 

i -i \ 

I 

1 

I 

—i i 

“ 2 

1 

-I 

i i 


\i 

-I 

-i -i / 


{Xi^ (r)X,v ( 0 ))\ 

{Yi, ( 0 )) 

{Yi, ( 0 )) 

(Xi^(r)r,v'( 0 )) J 


(71) 


Depending on the symmetries of the cluster Hamilto¬ 


nian, some terms in 0 may be zero at all time and can 
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be removed from the computation for speed-up or used 
to monitor possible errors coming from noise or other 
sources. 


possible that the method can be extended to simulate 
non-equilibrium processes [49] by measuring the Keldysh 
matrices and G^. 


CONCLUSION 


ACKNOWLEDGMENTS 


We have outlined a method to compute different ob¬ 
servables of the FHM using a quantum computer. It syn¬ 
thesizes and builds mainly on the work of [ansi uni Ea 
ETj. Provided that the lattice can be divided into clus¬ 
ters (with Lc spin-^ orbitals) which are coupled only with 
one-body hopping terms, section [TT| reviewed how a vari¬ 
ational principle for the grand canonical potential of the 
model can be used to approximate the self-energy of the 
lattice Hamiltonian and account for possible long-range 
ordering effects. A similar construction where a func¬ 
tional would also integrate an interaction across clusters 
[23] could also be considered. 

The formalism to define a cluster was reviewed in sec- 
tion |III| through the form of an example 2D lattice divided 
in 2 X 2 clusters for which a few order parameters like anti¬ 
ferromagnetism and superconductivity can be described 
and observable quantities computed. However, assuming 
no spin, spatial or electron-hole symmetries in the cluster, 
up to AL‘1 variational terms can be defined. The nature of 
the saddle-point problem that needs to be solved numer¬ 
ically is detailed and the bottleneck is shown to be the 
diagonalization and the simulation of the cluster which 
have to be solved for several variational parameters. 

The scaling and solution methods for a given cluster 


are detailed in section IV The memory scaling is known 
to be very bad on classical computers as the dimension 
of the Hilbert space of a cluster scales as 4^^ in the num¬ 
ber of orbitals. A method which assumes some way of 
creating a Gibbs state at low temperature on a quantum 
computer is presented. It is shown that there are 41/^ 
time correlation functions that need to be measured each 
round of the saddle-point optimization problem. The 
Bravyi-Kitaev transformation is known to significantly 
improve the scaling of classical algorithm in the case 
where the number of electrons is conserved by the Hamil¬ 
tonian [43| but a similar ansatz may also improve the 
method presented in this paper (by dividing the Hilbert 
space in even/odd occupation blocks for example). 

This algorithm provides a novel way to simulate com¬ 
plex materials at the electronic level and study new ques¬ 
tions without knowing the answer in advance. However 
some aspects could be improved. Notably, it is not fully 
clear whether the transformation on the Gibbs state be 
conditionally reversed after a measurement in such a way 
that the state can be reused. The back-action of the cor¬ 
relation function measurement may prevent the recycling 
of the Gibbs state. Also, it may be possible to estimate 
the errors of the algorithm by simulating a known sys¬ 
tem and comparing with analytical results (for example 
one could simulate the well-known tight-binding model 
to benchmark the quantum algorithm). Finally, it is 


We are grateful to David Senechal for the very helpful 
discussion. This work was supported by the European 
SCALEQIT program and Saarland University. 


Appendix A: Numerical example on the ID chain 

The simplest experimental implementation of the vari¬ 
ational procedure on a quantum computer would corre¬ 
spond to solving a simple ID tight-binding chain. With 
a minimum cluster of Lc = 2 sites (labeled “I” and “2”) 
each with 2 electrons (spin-up and spin-down), a 5-qubit 
quantum computer would be sufficient to extract the cor¬ 
relations functions ( [6^ . This section shows in detail how 
the formalism of subsection |IV C| can be used to compute 
the band structure and its occupation for the ID chain at 
arbitrary /i and T. The simulation was restricted only to 
a chemical variational potential ji' and a simple pairing 
potential A' which is expected to be zero in the case of 
one dimension. 


1. Finding the saddle-point of the self-energy 
functional 


First, the saddle point f I of equation (26) must 


be found. This is done through the following sequence: 


I. Choose a point 


hi[±h\ 

V J 

parameter). 


and 



and its neighbors 
(with h a small 


2. On a quantum computer, measure the retarded 

Nambu Greenes function (r, , A') of the clus- 

ter for the points of step 1 (as described in section 

3. Numerically compute the square of the gradient 
(26). If the modulus of the gradient is smaller 

f nN 

than some threshold cq, stop and assign I j ~ 
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4. Using a numerical Newton-Raphson method m, 
pick the next point ( ^ and loop over to step 

1 . 


Once the saddle-point is known, (r,A'^) is mea¬ 
sured and properties like the spectral density of the lat¬ 
tice can be approximated. 


2. Measuring and calculating the retarded Green’s 
function of the cluster 


The retarded Nambu Green’s function is measured on 
a discrete time domain = nAr where n is an inte¬ 
ger between 0 and rimax and Ar is a small time interval 
(^max = 2000 and Ar = 0.05 in this example) such that 
^max = ^maxAr. The matrix form of G'^ clearly shows 
that the number of correlation functions (c^ (r) cj, (0)) 
scales as AL‘i\ 


G {jn) = -iO (Tn) 


(. 

(cit (Tn)T (0)) 

(cit {Tn)cl^ (0)\ 

(Cit(T-n)cu (0)) 

(Cit iT-n)C2l (0)) 

g 


(C2t (Tn)4t W) 

(C2t {Tn)clt (0)) 

(C2t (^n)cU (0)) 

(c2t {Tn)c2l (0)) 



(A ijn)c\^ (0)) 

(A K)c2t (0)/ 

(A (^n)cu (0)) 

(T (Tn)c2l (0)) 


(' 

(^«) 4t (®)) 

( 44 . (Ai) C2t 

\4l (T'n)cu (0)) 

(T ('rn)C2i (0)) 

' ) 


It is then Fourier transformed on a discrete frequency 
domain ujm = mAuj between —cjmax and cjmax chosen 
such that C^max 


= 2 ^ and Aw = 2 ^ 

^ I max 


i /i? 


A '^max 


(Tn). 


(A2) 


n=0 


The numerical G'^ (cj) can then be used to compute the 
lattice-perturbed Green’s function Q ^k, (see equation 

([ 4 T])) and variou s properties of the latt ice as detailed in 
subsection III G The exact mapping of on the quan¬ 
tum computer is done through the Jordan-Wigner trans¬ 
formation 


■"It 


^2t 


= I(g)I( 
= I(g)I( 


) I G cr+ 


) < 


= I G cr+ 
t 


) az 


) az 


CTz 


) (Jz 


) CFz 


(Al) 


(A3) 


^2t — ' 

Using this transformation, all component of the Hamil¬ 


tonian 1-L' of the cluster (15) are mapped to a 4-qubit 
Hilbert space: 


'Wfh = -t {c\^C2^ + 4t^U + 4t^2t + “ U (ni^ni; + n2^n2i) 

= -t (I (g) I (g) (cr_ (g) (7+ -h cr+ (g) cr_) + (cr_ (g) (7+ -h (7+ (g) cr_) (g) I (g) I) 

— U (I (g) (7^ (g) I (g) (g) I (g) (g) I) , 

Rpair — A T T ^2^4^, T C2t^2t^ 

= A' (I (g) (cr+ 0 (7z 0 cr^ + (7- 0 (Tz <S) cr_) -f (( 7 + (g) cr;^ (g) cr+ + cr_ (g) cr;^ (g) cr_) (g) I) , 


(A4) 


(A5) 


^local = U' (^It + ^2t + + ^2t) 

= /i'(I(g)I(g)I(g)crn-hII(g)II(g)cr^(g)II + II(g)cr^(g)II(g)II + cr^(g)II(g)II(g)II). 


(A6) 


It can be noticed that the standard Fermi-Hubbard chemical potential can be implemented with single qubit 
term requires gates between two qubits, the variational 
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gates but the pairing terms need operations over several Finally the operators that are applied in the phase es- 
qubits to maintain the statistics of the fermions. The timation part of the algorithm and are required in the 
perturbation matrix (33) is given explicitly by reconstruction of (Al) are given by the following trans¬ 

formations: 

/ —/i + /i' e(k^-\-t —A' 0 \ 




/ e + t 
e* + t —fi + ii' 

-A' 0 


0 


V 


-A' 


.•(i) 


0 

-A' 

—e — t 

— t ji — ii' 

(A7) 


/ 


Xlt 

II 

+ 

II 

(g) I (g) I G (Ja,, 

Ut = 

—i 1 

(cit 

- T) 

= (g) I (g) I (g) 

X2t 

= C2t + 4i- = 

(g) I (g) (Ja; (g) az, 

F2t = 

—i 1 

(c2t 

- C2t) 

= ay 0 az, 


= Ci^ + cJj, = 

CFx^ (Jz^ 
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Figure 7. Measured probabilities for different and at 
different times. In this case the cluster parameters are Lc — 2^ 
t = 1, = A' = /i' = 0 and T = 0.1. 


The procedure highlighted in subsection |IIIB| is then 
followed to compute the CPT Green’s function and the 
desired properties of the system. 


3. Simple tight-binding model 

The tight-binding model U = 0 is investigated using 
the methods of this paper. The goal is to show that 
the method can accurately simulate well known simple 
models through the intermediate results it produces. 


Figure 8 . Non-zero correlation functions computed from the 
results of figure The function was regularized with a 
term to remove the fast oscillations of the Fourier transform 
arising from the finiteness of the time domain, 77 = ^ was 
used in this case. 


In figure the measured value of (Ad = 1,t) is 
shown for the simplest case of a 2-site tight-binding clus¬ 
ter. In this case the model generates simple oscillations 
as no decoherence is included. 

In figure the Green’s functions (r) computed 

from equation 0 are shown. Notice that the time- 
dependent Green’s functions were regularized with an 
decaying exponential in order to remove the fast os¬ 
cillations coming from the convolution of the frequency- 
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Figure 9. Real and imaginary parts of the frequency- 
dependent Green’s functions arising from the correlation func¬ 
tions measured in figure 



20 


10 

0 


-10 



Figure 10. Potthoff functional Q for different variational pa¬ 
rameters jji' and A' of a cluster of size Lc — 2 with parameters 
t=l,f/ = 0, = 0 and T = 1. The cross marks the saddle 


point at 




0.0046 

0 


Figure 11. Electron momentum-frequency distribution 
A{k,cj) for a lattice with parameters t=l, 1/ = 0, /i = 0 
and T = 1. The cluster used had Lc — 2 site and the saddle- 
point is the same as in figure The dashed line is at the 
chemical potential. 



Figure 12. Electron momentum distribution N {k) for differ¬ 
ent chemical potentials fi and temperatures T with = 0. 
The solid lines are the results from the numerical simulation 
of the quantum algorithm using time steps of size dr — 0.02 
up to Tmax = 200 whilc the dashed lines come from an imagi¬ 
nary frequency summation. 


dependent Green’s function with the sine ( ^ 
involved in finite time measurements. This regularizing 
term is not decoherence, but it could model a uniform 
depolarizing rate r] in the quantum processor. This rate 
would actually contribute to the width of the frequency- 
dependent Green’s function. 

In figure M the Fourier transformed {uj) are 

shown for the simple tight-binding cluster. Only two 
peaks are present and their width is determined by 77 
and the time domain used to measure the correlation 
functions. 

Figure shows an example of the Potthoff functional 
and its saddle point for a small ID cluster. 
As expected for this simple model, the saddle point is 


almost at the origin, the small deviation comes from the 
low finite temperature. At the saddle point, the average 
occupation of each state is (n) = 0.5 as is expected. At 
the saddle-point the spectral density of the full lattice 
can be computed. 

Figure pT] shows the spectral density A (/c, ci;)computed 
from equation ( [4^ for 50 clusters of size Lc = 2 in a 
simple tight binding model at relatively high temperature 
T = 1. The cosine band is fill above the Fermi level 
because of the high temperature. 

Figure[^shows that the simulation yields the expected 
physics of the tight-binding model at finite temperature. 
The ground state is indeed a ID Fermi sea in the elec- 
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k/TT 


Figure 13. Electron momentum-frequency distribution 
A (k,cj) for a lattice with parameters t=l,f/ = 4,/i = —2 
and T = 0.1 . The cluster used had Lc — 2 site and the 

. The dashed line is at 

the chemical potential. 


saddle-point is at 


A' 


-2 

0 


tronic momentum distribution (47) whose width is in¬ 
creased with the chemical potential and broadened by 
increased temperature. The loss of accuracy in the simu¬ 
lation is attributed to the sampling method and the accu¬ 
racy of the Fermi distribution on the discrete frequency 
domain computed from the measured time series. 

Finally figure 13 shows the spectral density A{k,uj) 
computed from equation for a cluster of size Lc = 2 in 
an attractive Hubbard chain 1/ = 4 at low temperature 
T = 0.1. The band is highly distorted by the interaction 
and the ground state is no longer a /c = 0 state. 

Extending these calculation for more complicated 
model is an easy task. A simple 2D model with a su¬ 
perconducting phase transition would require 4 sites and 
8 electrons, so a 9-qubit quantum computer would be re¬ 
quired to measure (r) in this case. It appears that 
the number of time points that need to be measured may 
become an issue as the systems become more complex. 
It would be interesting to know if there exist sampling 
methods as efficient as imaginary frequency summation 
methods where only ^ 100 points need to be mea¬ 
sured in order to achieve a high numerical accuracy in 
the computation of the Fermi function even for compli¬ 
cated electronic structures. For example, a cost function 
over several models could be used to extract the Green’s 
function using fewer measurements. Alternatively, mea¬ 
suring forward finite difference time derivatives close to 
r = 0 to get the coefficients of the moment expansion 
of equation (41) could also work. Indeed, the correlation 
functions (64) can be rewritten as 


R - 




Pi P2 Ps 


1 




1 


Pqc 


Figure 14. Detailed circuit to prepare an approximate Gibbs 
state Pqc ~ pcibbs following m- The simulated inverse tem- 
perat ure /3 is related to the measurement of 5=,= by equation 
( B1Q| . The initial state of R and Q is taken to be the zero 
state then the Hadamard gate is applied 

on R and Q is transformed (non-unitarily) to the fully mixed 
state Thcu q controlled-!/ operations are ap- 

_ . _Hq_ 

plied, where the notation Ut — and U — e ^H^olloo with 
Ro — R' Rb- An inverse quantum Fourier transform is 
applied on register R and the string s* is read from the first 
q qubits. Register S is then left in a simulated Gibbs state 
Pqc- 


where the moments are given by 








= (At)-* EEo (-!)’■ ( " ) ((" - P 

rO (At) 


(AlO) 

which could be approximated experimentally by forward 
finite differences (higher order finite differences could also 
be used). 


Appendix B: Preparation of a Gibbs state 

A digital method to prepare Gibbs states in a quan¬ 
tum computer is reviewed and shown adequate for a vari¬ 
ational solver. The goal is the make this document self- 
contained in the sense that action of the quantum com¬ 
puter can be fully defined. 

Here is the summary of the method, as given in m , to 
prepare the Gibbs state required to simulate the correla¬ 
tion function of the cluster. In addition to the simulated 
system Hamiltonian R'^ a bath Hamiltonian Rb is re¬ 
quired such that the total uncoupled system is 


5=0 


s) 

flU 


(A9) 


Ho = H'+Hb (B1) 


with eigenvalues 



and energy eigenvectors 
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E 


( 0 ) 


^ I. The bath (first part of the register Q in figure 

pI is assumed to be a collection of m uncoupled spin-^ 
with energy splitting 77 : 


= (B2) 

i=i 

A small interaction V is allowed such that the total cou¬ 
pled system Hamiltonian is 


ntot = Ho + V (B3) 


with eigenvalues {Ej^} and energy eigenvectors {\E]^)}. 
The procedure is the following (see figure [T^ 


1. Initialization, r Hadamard gates H are applied on 
the qubits of register R and the register Q is relaxed 
in the fully mixed state of (B3) sueh that 


||Htot Iloo inverse temperature P is de¬ 

termined by E and Arect- The final state in the 
register Q is now 

Pqc = TrRps 


oc Eti (EtTiir • 

(B 8 ) 

One of the rectangular states contained in the ini¬ 
tial fully mixed state is selected upon measurement. 
For appropriately chosen parameters, the state in 
register S is approximately a Gibbs state of the clus¬ 
ter Hamiltonian. 


The algorithm outputs a reduced state = TtbPqc 


/^Gibbs ~ chanucl S, where /3 = ^ is the 

inverse temperature. Assuming a bath of the form (B2) 

with energy scale p = ||H'||^, the really implies 

the following condition 





s,s'=0 k=l 


(B4) 


^ {pQC^ Pcibbs) ^ ^ 


ln( 2 "-«)^ eX+g|l'«'lloo + - 


2r-q-2 


where d = is the total dimension of the 

system plus bath. This is equivalent to preparing 
the coupled system + bath at infinite temperature. 


2. Partial quantum phase estimation, r controlled-U 
operation are followed by an inverse Fourier trans- 

i ^0 

form on R. Note that U = e ^'^oiioo ^ where Ho is 
the uncoupled Hamiltonian (Bl). After this phase 
estimation part, the state in the computer is 


^ 2^-1 d 

P 2 = ^ ^ ^ ^ ^ (pk) (pk) |«5) {s I 0 \Tk) {Tk\ 

s,s'=0 k=l 

(B5) 

where pk = fuFRl — 


I I — ^27Ti{2^ip-s) 

(t) = ^ 2 — g27ri((/?-2-^s) 

The controlled evolution of the full system dephases 
different distributions of eigenvalues contained in 
the fully mixed state. 

3. Measurement. The first q qubits of R are measured. 
A binary string 5 ^ (length q) is obtained 


-|-l)Z\rect=t= d 

P3 OC E E“* {<Pk)a:, i^k)\s) {s'\®\Ek){Ek\ 

S ,s'=S ^/\rect* k = l. 

(B7) 

where Arect* = 2 ^^ ^ is the number of states of the 
ancillary register R compatible with the measure¬ 
ment. The width of the rectangular state that is pre¬ 
pared is determined by Arect = ||Htot ||00 2“"^Arect*- 
The energy of the rectangular state is E = 


+ 1 + C 

(B9) 

where T>(-,-) is the trace distance and C is a constant 
exponentially small in m. The effective inverse temper¬ 
ature is in the interval [^ — 6/3 F 6/3] with 




4 

p 


- 2-^5* 



PbIl;T 


(BIO) 


Since s^ G [ 0 , 2 *^ — 1 ], the inverse temperature of the gen¬ 
erated Gibbs state can reach negative values in princi¬ 
ple (physically corresponding to a state with an inverted 
population). The uncertainty on the temperature of the 
Gibbs state is bounded by 


-5/3 < 4G(i + 


ll«'IL 

II«b||, 


(Bll) 


_ o2—q /A 

“ ^ y m\\n 


E ( 


1 + 


4 . 


VmX ) 

At least q qubits are needed according to the rule 


q> 


-l0g2 



Spr] \ 

ll«'IU 
IIBbIU / 


+ 2 


(B12) 


and the average number of runs required to achieve some 
inverse temperature is 


ttruns < 2 ^ 




(B13) 


This last bound is a worst-case scenario as finding 
the ground state of the Fermi-Hubbard is in general a 
QMA — hard problem. 
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